# Title: Public Trust and Collaborative Governance: An Instrumental Variable Approach
# Yixin Liu 2022
# Public Management Review
# PMR Replication file for Study 1

library(ggplot2) # plot
library(ggridges) # overlapping line plots
library(pwr) # power
library(ggpubr) # plot
library(effsize) # effect size
library(multiwayvcov) # cluster se
library(lmtest) # cluster se
library(egg) # ggplot(theme_article)
library(cowplot) # ggplot output
library(Rmisc) # plot stat table function
library(stargazer) # latex
library(xtable) # latex
library(texreg) # latex
library(ivreg)
library(standardize)
library(Rmisc) # summaryse function

####Function####
element_textbox <- function(...) {
  el <- element_text(...)
  class(el) <- c("element_textbox", class(el))
  el
}

element_grob.element_textbox <- function(element, ...) {
  text_grob <- NextMethod()
  rect_grob <- element_grob(calc_element("strip.background", theme_bw()))
  
  ggplot2:::absoluteGrob(
    grid::gList(
      element_grob(calc_element("strip.background", theme_bw())),
      text_grob
    ),
    height = grid::grobHeight(text_grob), 
    width = grid::unit(1, "npc")
  )
}

# Setup working directory
WD = getwd()
setwd(WD)

####Data####
exp<-read.csv("study1_data.csv")

####Trust####
# ANOVA test
trust.aov <- aov(trust~rand, data = exp)
summary(trust.aov)

####Figure 2####
# left panel
trust.df<-summarySE(exp, measurevar="trust", 
                    groupvars="rand",
                    na.rm = T)

trust.df$rand<-as.factor(trust.df$rand)
trust.df$rand<-ordered(trust.df$rand,level=c("Corrupt","Control","Honest"))

trust.plot<-ggplot(trust.df,aes(x=rand, y=trust, group=1)) +
  geom_line(linetype = "dashed") +
  geom_point(aes(size=3)) +
  geom_pointrange(aes(ymin=trust-ci, ymax=trust+ci)) +
  scale_y_continuous(lim=c(30,80)) +
  labs(x = element_blank(),
       y = element_blank(),
       title = "Trust in Government")+
  theme_article()

trust.plot<-trust.plot + theme(legend.position = "none",
                               axis.title.y = element_text(size = 12),
                               axis.text.x=element_text(size = 12),
                               axis.text.y=element_text(size = 10),
                               strip.background = element_blank(),
                               strip.text.x = element_blank(),
                               axis.ticks.x = element_blank(),
                               plot.title = element_textbox(
                                 hjust = 0.5, margin = margin(t = 5, b = 5)))

trust.plot

# right panel
private.df<-summarySE(exp, measurevar="private", 
                      groupvars=c("rand"),
                      conf.interval = 0.95,
                      na.rm = T)

community.df<-summarySE(exp, measurevar="community", 
                        groupvars=c("rand"),
                        conf.interval = 0.95,
                        na.rm = T)

coproduce.df<-summarySE(exp, measurevar="coproduce", 
                        groupvars=c("rand"),
                        conf.interval = 0.95,
                        na.rm = T)

names(private.df)[names(private.df) == "private"] <- "dv"
private.df$group<-"Public-Private Partnerships"

names(community.df)[names(community.df) == "community"] <- "dv"
community.df$group<-"Public-Citizen Partnerships"

names(coproduce.df)[names(coproduce.df) == "coproduce"] <- "dv"
coproduce.df$group<-"Willingness to Coproduce"

df.dv<-rbind(private.df,community.df,coproduce.df)

df.dv$group <- factor(df.dv$group, 
                      levels = c("Public-Citizen Partnerships", 
                                 "Willingness to Coproduce",
                                 "Public-Private Partnerships"))

dv.plot<-ggplot(df.dv,aes(x=rand, y=dv, group=group, shape=group)) +
  geom_line(linetype = "dashed") +
  geom_point(aes(shape=group), size = 4) +
  geom_pointrange(aes(ymin=dv-ci, ymax=dv+ci)) +
  scale_y_continuous(lim=c(50,80)) +
  labs(x = element_blank(),
       y = element_blank(),
       title = "Percieved Collaborative Governance")+
  theme_article()

dv.plot<-dv.plot + theme(legend.position = c(0.75, 0.15),
                         legend.direction = "vertical",
                         legend.title = element_blank(),
                         legend.background = element_blank(),
                         legend.text = element_text(size = 12),
                         strip.background = element_blank(),
                         strip.text.x = element_blank(),
                         axis.ticks.x = element_blank(),
                         axis.title.y = element_text(size = 12),
                         axis.text.x=element_text(size = 12),
                         axis.text.y=element_text(size = 10),
                         plot.title = element_textbox(
                           hjust = 0.5, margin = margin(t = 5, b = 5)))

dv.plot

# left + right panel
combine.plot<-ggarrange(trust.plot,
                        dv.plot,
                        ncol = 2, nrow = 1)
combine.plot<-plot_grid(combine.plot)
combine.plot
#ggsave(file = "figure2.jpg", width = 12, height = 5, dpi = 666)

####Table 2####
# stage 1 (model 1)
trust1<-lm(trust~trustprime,data=exp)
summary(trust1)

# stage 2 (model 2)
iv1 = ivreg(private~trust
            | trustprime, data = exp)
summary(iv1)

# stage 2 (model 3)
iv2 = ivreg(community~trust 
            | trustprime, data = exp)
summary(iv2)

# stage 2 (model 4)
iv3 = ivreg(coproduce~trust 
            | trustprime, data = exp)
summary(iv3)

# regression table
main1<-list(trust1,iv1,iv2,iv3)
texreg(main1,
       custom.model.names=c("(1)","(2)","(3)","(4)"),
       custom.coef.map=list("trust"="Trust in Government",
                            "trustprime"="Government Integrity",
                            "(Intercept)"="Constant"),
       custom.note = paste("\\item %stars. Government integrity is measured 
                           by coding experimental conditions as: 
                           1 = Honest; 0.5 = Control; 0 = Corrupt."),
       include.adjrs = FALSE,include.nobs = TRUE, single.row = FALSE,
       fontsize="footnotesize", caption.above=TRUE,digits=2,dcolumn=TRUE,
       caption="2SLS Estimate in Study 1")

####Table 3####
# subgroup separate instrument
exp.honest<-subset(exp, exp$trustprime!=0)
exp.honest$trustprime<-ifelse(exp.honest$trustprime==1,0,1)
exp.corrupt<-subset(exp, exp$trustprime!=2)
exp.corrupt$trustprime<-ifelse(exp.corrupt$trustprime==1,0,1)

# Panel A
# stage 1 (model A1)
honest1<-lm(trust~trustprime,data=exp.honest)
summary(honest1)

# stage 2 (model A2)
iv1.honest = ivreg(private~trust
                   | trustprime, data = exp.honest)
summary(iv1.honest)

# stage 2 (model A3)
iv2.honest = ivreg(community~trust
                   | trustprime, data = exp.honest)
summary(iv2.honest)

# stage 3 (model A4)
iv3.honest = ivreg(coproduce~trust
                   | trustprime, data = exp.honest)
summary(iv3.honest)

# Panel B
# stage 1 (model B1)
corrupt1<-lm(trust~trustprime,data=exp.corrupt)
summary(corrupt1)

# stage 2 (model A2)
iv1.corrupt = ivreg(private~trust
                   | trustprime, data = exp.corrupt)
summary(iv1.corrupt)

# stage 2 (model A3)
iv2.corrupt = ivreg(community~trust
                   | trustprime, data = exp.corrupt)
summary(iv2.corrupt)

# stage 3 (model A4)
iv3.corrupt = ivreg(coproduce~trust
                   | trustprime, data = exp.corrupt)
summary(iv3.corrupt)

# regression table
# Panel A
honest1<-list(honest1,iv1.honest,iv2.honest,iv3.honest)
texreg(honest1,
       custom.model.names=c("(1)","(2)","(3)","(4)"),
       custom.coef.map=list("trust"="Trust in Government",
                            "trustprime"="Honest",
                            "(Intercept)"="Constant"),
       custom.note = paste("\\item %stars. Honest is measured by 
                            coding experimental conditions as: 1 = Honest; 0 = Control."),
       include.adjrs = FALSE,include.nobs = TRUE, single.row = FALSE,
       caption.above=TRUE,digits=2,dcolumn=TRUE,
       caption="2SLS Estimate Using Honest and Corrupt as Separate Instruments in Study 1 (Panel A)")

# Panel B
corrupt1<-list(corrupt1,iv1.corrupt,iv2.corrupt,iv3.corrupt)
texreg(corrupt1,
       custom.model.names=c("(1)","(2)","(3)","(4)"),
       custom.coef.map=list("trust"="Trust in Government",
                            "trustprime"="Corrupt",
                            "(Intercept)"="Constant"),
       custom.note = paste("\\item %stars. Corrupt is measured by 
                            coding experimental conditions as: 1 = Corrupt; 0 = Control."),
       include.adjrs = FALSE,include.nobs = TRUE, single.row = FALSE,
       caption.above=TRUE,digits=2,dcolumn=TRUE,
       caption="2SLS Estimate Using Honest and Corrupt as Separate Instruments in Study 1 (Panel B)")

####Appendix A####
#####Appendix A.2####
# study 1 sample
vars<-data.frame(exp$sex,
                 exp$white,exp$black,exp$hispanic,exp$asian,exp$other,
                 exp$ideology,
                 exp$age,
                 exp$education,
                 exp$income,
                 exp$gw,exp$gwh,exp$climatebelief)

Mean=sapply(vars, mean, na.rm = TRUE)
SD=sapply(vars, sd, na.rm = TRUE)
Min=sapply(vars, min, na.rm = TRUE)
Max=sapply(vars, max, na.rm = TRUE)

female.aov <- aov(sex ~ rand, data = exp)
summary(female.aov)
female<-summary(female.aov)[[1]][,5][1]

wht.aov <- aov(white ~ rand, data = exp)
summary(wht.aov)
wht<-summary(wht.aov)[[1]][,5][1]
blk.aov <- aov(black ~ rand, data = exp)
summary(blk.aov)
blk<-summary(blk.aov)[[1]][,5][1]
his.aov <- aov(hispanic ~ rand, data = exp)
summary(his.aov)
his<-summary(his.aov)[[1]][,5][1]
asian.aov <- aov(asian ~ rand, data = exp)
summary(asian.aov)
asian<-summary(asian.aov)[[1]][,5][1]
other.aov <- aov(other ~ rand, data = exp)
summary(other.aov)
other<-summary(other.aov)[[1]][,5][1]

ideo.aov <- aov(ideology ~ rand, data = exp)
summary(ideo.aov)
ideo<-summary(ideo.aov)[[1]][,5][1]

age.aov <- aov(age ~ rand, data = exp)
summary(age.aov)
age<-summary(age.aov)[[1]][,5][1]

educ.aov <- aov(education ~ rand, data = exp)
summary(educ.aov)
educ<-summary(educ.aov)[[1]][,5][1]

inc.aov <- aov(income ~ rand, data = exp)
summary(inc.aov)
inc<-summary(inc.aov)[[1]][,5][1]

gw.aov <- aov(gw ~ rand, data = exp)
summary(gw.aov)
gw<-summary(gw.aov)[[1]][,5][1]

gwh.aov <- aov(gwh ~ rand, data = exp)
summary(gwh.aov)
gwh<-summary(gwh.aov)[[1]][,5][1]

climate.aov <- aov(climatebelief ~ rand, data = exp)
summary(climate.aov)
climate<-summary(climate.aov)[[1]][,5][1]

# sample table
des_stat<-matrix(NA,length(Mean),5)
colnames(des_stat)<-c("Mean","SD","Min","Max","Randomization Check (P-value)")
rownames(des_stat)<-c("Female","White","Black","Hispanic","Asian","Other",
                      "Ideology","Age","Education","Income", "GW", 
                      "GWH", "Climate Change Beliefs")

des_stat[,1]<-Mean
des_stat[,2]<-SD
des_stat[,3]<-Min
des_stat[,4]<-Max
des_stat[,5]<-c(female,wht,blk,his,asian,other,ideo,age,educ,inc,gw,gwh,climate)

des_stat<-round(des_stat,digits = 2)

# write table
xtable(des_stat, digits=2)

#####Appendix A.3####
# Placebo test
iv.p1 = ivreg(placebo1~trust 
              | trustprime, data = exp)
summary(iv.p1)

iv.p2 = ivreg(placebo2~trust 
              | trustprime, data = exp)
summary(iv.p2)

placebo<-list(iv.p1,iv.p2)
texreg(placebo,
       custom.model.names=c("(1)","(2)"),
       custom.coef.map=list("trust"="Trust in Government",
                            "(Intercept)"="Constant"),
       custom.note = paste("\\item %stars. "),
       include.adjrs = FALSE,include.nobs = TRUE, single.row = FALSE,
       fontsize="footnotesize", caption.above=TRUE,digits=2,dcolumn=TRUE)

#####Appendix A.4####
trustcat<-lm(trust~rand,data=exp)
summary(trustcat)

cat1<-lm(private~rand, data = exp)
summary(cat1)

cat2<-lm(community~rand, data = exp)
summary(cat2)

cat3<-lm(coproduce~rand, data = exp)
summary(cat3)

cat<-list(trustcat,cat1,cat2,cat3)

# regression table
texreg(cat,
       custom.model.names=c("(1)","(2)","(3)","(4)"),
       custom.coef.map=list("randHonest"="Honest",
                            "randCorrupt"="Corrupt",
                            "(Intercept)"="Constant"),
       custom.note = paste("\\item %stars. "),
       include.adjrs = FALSE,include.nobs = TRUE, single.row = FALSE,
       caption.above=TRUE,digits=2,dcolumn=TRUE)


